# -*- coding: utf-8 -*-
"""
海洋涡旋三维分层提取与速度场融合脚本
====================================
功能概述：
    基于DGGS（离散全球网格系统）六边形网格，整合涡旋原始数据与三维海洋速度场数据，
    实现涡旋的表层初始化、分层外扩提取、旋转特性验证，并输出二维（表层）与三维（全深度）涡旋结果。
    核心目标是过滤虚假涡旋（深度层级<3），保留具有物理意义的三维涡旋结构。

依赖库：
    - numpy: 数值计算与数组处理
    - netCDF4: NetCDF文件读写
    - math: 基础数学运算（如三角函数、哈维正弦公式）
    - tqdm: 处理进度可视化
    - os: 文件路径与目录操作
    - re: 正则表达式（提取文件名日期标识）
    - PyDGGS: 离散全球网格系统工具（六边形编码/解码、邻域查询）

全局常量定义
------------
- Rossby_Radius: 罗斯贝半径阈值（单位：km），用于过滤过小涡旋（<该值的涡旋不处理）
- dggs_level: DGGS网格分辨率等级（等级越高，网格单元越小）
- grid_type: DGGS网格类型（ISEA4H六边形网格）
- element_type: DGGS网格元素类型（单元格）
- topo_type: DGGS拓扑查询类型（相对单元格查询，用于获取邻域）
- kernel: DGGS内核对象，用于坐标编码/解码、邻域查询等核心操作
"""

import numpy as np
import netCDF4 as nc
import math
from tqdm import tqdm
import os
import re
import PyDGGS

# 全局常量配置（根据业务需求调整）
Rossby_Radius = 10  # 罗斯贝半径阈值（单位：km），过滤过小涡旋
dggs_level = 8  # DGGS网格分辨率等级
grid_type = PyDGGS.DGGSGridType_DGGS_ISEA4H  # DGGS网格类型：ISEA4H六边形
element_type = PyDGGS.DGGSElementType_Cell  # DGGS元素类型：单元格
topo_type = PyDGGS.DGGSTopoType_RelativeCells  # DGGS拓扑查询类型：相对单元格
kernel = PyDGGS.CreateKernelObj(grid_type)  # 初始化DGGS内核对象


def get_hash_key(code):
    """
    将DGGS单元格编码转换为哈希键（元组形式），用于字典索引

    Args:
        code (PyDGGS.DGGSElementCode): DGGS单元格编码对象，包含row（行）、col（列）、basePartition（基础分区）属性

    Returns:
        tuple: 哈希键，格式为 (code.row, code.col, code.basePartition)
    """
    return (code.row, code.col, code.basePartition)


class EddyHybridIntegration:
    """
    涡旋-速度场混合整合类

    核心职责：
        1. 初始化DGGS六边形网格单元格数据
        2. 处理涡旋原始数据，完成表层（depth_level=0）涡旋ID分配
        3. 融合三维速度场数据，实现涡旋的分层外扩与旋转特性验证
        4. 计算各深度层级涡旋轮廓、中心与半径
        5. 过滤虚假涡旋（深度层级<3），保存二维（表层）与三维（全深度）结果

    属性说明：
        - params (dict): 涡旋验证参数（如速度连续性、旋转角度容忍阈值等）
        - depth_levels (int): 深度层级总数（固定30层，索引0~29，与速度场深度维度一一对应）
        - hex_resolution (int): DGGS网格分辨率等级（与全局dggs_level一致）
        - grid (dict): 单元格数据字典，键为 (hex_code, depth_level)，值为单元格详细信息（经纬度、速度、涡旋ID等）
        - eddy_cells (dict): 涡旋单元格集合，键为eddy_id，值为涡旋各深度层级信息（中心、单元格、轮廓、半径）
        - false_eddies (list): 虚假涡旋ID列表（深度层级总数<3的涡旋）
        - eddy_data (netCDF4.Dataset): 涡旋原始数据NetCDF对象
        - velocity_data (netCDF4.Dataset): 速度场数据NetCDF对象
        - is_cyclonic (bool): 涡旋类型标识（True=气旋，False=反气旋）
        - failure_stats (dict): 涡旋验证失败统计（按失败原因分类：数据质量、速度大小、旋转方向等）
        - eddy_params (dict): 涡旋原始元数据（从NetCDF文件读取的变量，如纬度、经度、半径等）
    """

    def __init__(self, is_cyclonic=True, hex_resolution=dggs_level):
        """
        类初始化函数

        Args:
            is_cyclonic (bool, optional): 涡旋类型（True=气旋，False=反气旋），默认True
            hex_resolution (int, optional): DGGS网格分辨率等级，默认与全局dggs_level一致
        """
        # 涡旋验证参数配置（可根据数据特性调整）
        self.params = {
            'Sv': 10,  # 速度连续性基础阈值
            'Sa': 180,  # 最大允许角度差基础阈值（单位：度）
            'Sae': 45,  # 旋转方向容忍阈值基础阈值（单位：度）
            'San': 7,  # 邻域验证数量阈值
            'Sd': 45,  # 深度连续性阈值（单位：度）
            'Sy': 180  # 对称性验证阈值（单位：度）
        }

        # 深度层级配置：固定30层，索引0~29，与速度场深度维度完全一一对应（无额外映射）
        self.depth_levels = 30
        self.hex_resolution = hex_resolution

        # 核心数据结构初始化
        self.grid = {}  # 键：(hex_code, depth_level)，值：单元格数据（经纬度、速度、涡旋ID等）
        self.eddy_cells = {}  # 键：eddy_id，值：涡旋各深度层级信息
        self.false_eddies = []  # 虚假涡旋ID列表（深度层级总数<3）
        self.eddy_data = None  # 涡旋原始数据NetCDF对象
        self.velocity_data = None  # 速度场数据NetCDF对象
        self.is_cyclonic = is_cyclonic  # 涡旋类型（气旋/反气旋）

        # 验证失败统计初始化（按失败原因分类）
        self.failure_stats = {
            'data_quality': 0,  # 数据质量不足（有效速度点<4）
            'velocity_magnitude': 0,  # 速度大小连续性不满足
            'rotation_direction': 0,  # 旋转方向与预期不符（气旋逆时针/反气旋顺时针）
            'velocity_structure': 0,  # 速度结构不连续（相邻点速度比过大）
            'symmetry': 0  # 旋转对称性/方向分布均匀性不满足
        }

    def _init_hex_cell(self, hex_code, depth_level, lat, lon):
        """
        初始化DGGS六边形单元格数据（若单元格未存在于grid中）

        Args:
            hex_code (tuple): DGGS单元格哈希键，格式 (row, col, basePartition)
            depth_level (int): 深度层级索引（0~29，与速度场深度维度一一对应）
            lat (float): 单元格纬度（单位：度）
            lon (float): 单元格经度（单位：度）
        """
        # 仅初始化未存在的单元格，避免覆盖已有数据
        if (hex_code, depth_level) not in self.grid:
            self.grid[(hex_code, depth_level)] = {
                'hex_code': hex_code,  # 单元格哈希键
                'depth_level': depth_level,  # 深度层级（与速度场深度维度对应）
                'eddy_id': -1,  # 涡旋ID（-1表示未分配）
                'u': np.nan,  # 东向速度（初始为NaN，后续从速度场填充）
                'v': np.nan,  # 北向速度（初始为NaN，后续从速度场填充）
                'is_cyclonic': self.is_cyclonic,  # 涡旋类型（继承类的涡旋类型）
                'lat': lat,  # 单元格纬度
                'lon': lon  # 单元格经度
            }

    def _get_hex_code(self, lat, lon):
        """
        将经纬度坐标转换为DGGS六边形单元格哈希键

        Args:
            lat (float): 纬度（单位：度）
            lon (float): 经度（单位：度）

        Returns:
            tuple: DGGS单元格哈希键，格式 (row, col, basePartition)
        """
        # 初始化DGGS地理坐标对象并赋值
        coord = PyDGGS.GeoCoord()
        coord.lon = np.float64(lon)  # 经度转换为float64（PyDGGS要求）
        coord.lat = np.float64(lat)  # 纬度转换为float64（PyDGGS要求）

        # 初始化DGGS单元格编码对象并配置参数
        cell_code = PyDGGS.DGGSElementCode()
        cell_code.level = self.hex_resolution  # 网格分辨率等级
        cell_code.elementType = element_type  # 元素类型（单元格）

        # 经纬度编码为DGGS单元格
        kernel.encode(coord, cell_code)

        # 转换为哈希键返回
        return get_hash_key(cell_code)

    def _get_neighbors(self, hex_code, radius=1):
        """
        获取DGGS六边形单元格的邻域（指定半径层级的相邻单元格）

        逻辑：
            1. 从输入hex_code重建DGGS编码对象
            2. 按半径层级迭代查询邻域（半径1=直接相邻，半径2=相邻的相邻，以此类推）
            3. 去重并返回邻域单元格的哈希键列表

        Args:
            hex_code (tuple): 目标单元格哈希键，格式 (row, col, basePartition)
            radius (int, optional): 邻域半径（层级），默认1（直接相邻）

        Returns:
            list: 邻域单元格哈希键列表（无重复）
        """
        # 半径<1时无邻域，返回空列表
        if radius < 1:
            return []

        # 从哈希键重建DGGS单元格编码对象
        code = PyDGGS.DGGSElementCode()
        code.row, code.col, code.basePartition = hex_code
        code.level = self.hex_resolution
        code.elementType = element_type

        all_neighbors = []  # 存储所有邻域单元格哈希键（去重）
        current_layer = [code]  # 当前层级的单元格（初始为目标单元格）

        # 按半径迭代查询邻域
        for _ in range(radius):
            next_layer = []  # 下一层级的单元格（当前层的邻域）
            for current_code in current_layer:
                # 初始化拓扑查询结果对象
                query_result = PyDGGS.QueryResult()
                # 查询当前单元格的相对邻域
                kernel.query(topo_type, current_code, query_result)

                # 遍历查询结果，提取邻域单元格哈希键
                for i in range(query_result.getSize()):
                    nbr_code = query_result.getGridCodeByIndex(i)  # 邻域单元格编码
                    nbr_hash = get_hash_key(nbr_code)  # 转换为哈希键

                    # 去重：仅添加未出现过的邻域
                    if nbr_hash not in all_neighbors:
                        all_neighbors.append(nbr_hash)
                        next_layer.append(nbr_code)  # 下一层级的查询基础

            current_layer = next_layer  # 更新当前层级为下一层级

        return all_neighbors

    def _calculate_contour(self, hex_codes, depth_level, eddy_id):
        """
        计算指定深度层级下，涡旋单元格集合的轮廓、中心与最大半径

        核心逻辑：
            1. 解码所有单元格的经纬度坐标
            2. 计算涡旋中心（经纬度平均值）
            3. 计算每个单元格到中心的地表距离（哈维正弦公式），确定最大半径
            4. 按16个扇形划分区域，取每个扇形内最远点作为轮廓点，保证轮廓平滑
            5. 按角度排序轮廓点，转换为DGGS哈希键返回

        Args:
            hex_codes (list): 目标深度层级的涡旋单元格哈希键列表
            depth_level (int): 深度层级索引（0~29）
            eddy_id (int): 涡旋ID（用于日志定位，无实际计算作用）

        Returns:
            tuple: (contour_hex_codes, center, max_radius)
                - contour_hex_codes (list): 轮廓单元格哈希键列表
                - center (tuple): 涡旋中心（lon, lat），若无单元格则为None
                - max_radius (float): 涡旋最大半径（单位：米），若无单元格则为0.0
        """
        # 无单元格时返回空结果
        if not hex_codes:
            return [], None, 0.0

        # 步骤1：解码所有单元格的经纬度坐标
        points = []  # 存储 (lon, lat) 坐标对
        for hex_code in hex_codes:
            # 重建DGGS单元格编码对象
            code = PyDGGS.DGGSElementCode()
            code.row, code.col, code.basePartition = hex_code
            code.level = self.hex_resolution
            code.elementType = element_type

            # 解码为地理坐标
            coord = PyDGGS.GeoCoord()
            kernel.decode(code, coord)
            points.append((coord.lon, coord.lat))

        # 步骤2：计算涡旋中心（经纬度平均值）
        lons, lats = zip(*points)  # 拆分经度、纬度列表
        center_lon = np.mean(lons)  # 中心经度（平均值）
        center_lat = np.mean(lats)  # 中心纬度（平均值）
        center = (center_lon, center_lat)

        # 步骤3：计算每个单元格到中心的地表距离（哈维正弦公式），确定最大半径
        distances = [
            self._haversine_distance(center_lat, center_lon, lat, lon)
            for lon, lat in points
        ]
        max_radius = max(distances) if distances else 0.0  # 最大半径（单位：米）

        # 步骤4：按16个扇形取最远点作为轮廓点（保证轮廓平滑）
        num_sectors = 16  # 扇形数量（越多轮廓越平滑，16为平衡值）
        sector_angles = np.linspace(-np.pi, np.pi, num_sectors + 1)  # 扇形角度范围（弧度）
        contour_points = []  # 存储轮廓点（lon, lat）

        for i in range(num_sectors):
            # 当前扇形的角度范围（左闭右开）
            sector_start = sector_angles[i]
            sector_end = sector_angles[i + 1]

            # 计算每个点到中心的方位角（弧度）
            point_angles = np.array([
                math.atan2(lon - center_lon, lat - center_lat)
                for lon, lat in points
            ])

            # 筛选当前扇形内的点
            sector_mask = (point_angles >= sector_start) & (point_angles < sector_end)
            sector_points = np.array(points)[sector_mask]  # 扇形内的点
            sector_distances = np.array(distances)[sector_mask]  # 扇形内点的距离

            # 取扇形内最远点作为轮廓点
            if len(sector_points) > 0:
                max_dist_idx = np.argmax(sector_distances)
                contour_points.append(sector_points[max_dist_idx])

        # 步骤5：按方位角排序轮廓点（保证顺时针/逆时针顺序）
        if contour_points:
            # 计算每个轮廓点到中心的方位角
            contour_angles = np.array([
                math.atan2(lon - center_lon, lat - center_lat)
                for lon, lat in contour_points
            ])
            # 按方位角排序（从小到大，保证顺序一致）
            sorted_indices = np.argsort(contour_angles)
            contour_points = [contour_points[i] for i in sorted_indices]

        # 转换轮廓点为DGGS单元格哈希键
        contour_hex_codes = [self._get_hex_code(lat, lon) for lon, lat in contour_points]

        return contour_hex_codes, center, max_radius

    def assign_eddy_ids(self, eddy_file, is_cyclonic):
        """
        处理涡旋原始NetCDF文件，完成表层（depth_level=0）涡旋ID分配与单元格初始化

        核心逻辑：
            1. 读取涡旋原始数据（经纬度、边界、半径、类型等）
            2. 过滤过小涡旋（有效半径<Rossby_Radius）
            3. 初始化涡旋中心、边界及边界邻域的表层（depth_level=0）单元格
            4. 计算表层涡旋轮廓、中心与半径，保存涡旋基础信息

        Args:
            eddy_file (str): 涡旋原始数据NetCDF文件路径
            is_cyclonic (bool): 涡旋类型（True=气旋，False=反气旋），覆盖类初始化的类型
        """
        # 更新涡旋类型（以输入参数为准）
        self.is_cyclonic = is_cyclonic

        # 读取涡旋原始NetCDF文件
        with nc.Dataset(eddy_file, 'r') as f:
            self.eddy_data = f  # 保存涡旋数据对象
            num_eddies = len(f.dimensions['obs'])  # 涡旋总数（obs为涡旋维度名）

            # 读取核心变量（经纬度、边界经纬度）
            eddy_lat = f['latitude'][:]  # 涡旋中心纬度（每个涡旋1个值）
            eddy_lon = f['longitude'][:]  # 涡旋中心经度（每个涡旋1个值）
            boundary_lat = f['effective_contour_latitude'][:]  # 涡旋边界纬度（每个涡旋NbSample个点）
            boundary_lon = f['effective_contour_longitude'][:]  # 涡旋边界经度（每个涡旋NbSample个点）
            boundary_sample_count = len(f.dimensions['NbSample'])  # 每个涡旋的边界采样点数

            # 读取涡旋类型（气旋/反气旋）：优先从cyc变量读取，无则按输入类型赋值
            try:
                eddy_type = f['cyc'][:]  # cyc变量：1=气旋，-1=反气旋
            except KeyError:
                # 无cyc变量时，按输入is_cyclonic统一赋值
                eddy_type = np.ones(num_eddies) * (1 if self.is_cyclonic else -1)

            # 保存涡旋原始元数据（所有含obs维度的变量）
            self.eddy_params = {}
            for var_name in f.variables:
                var = f[var_name]
                if 'obs' in var.dimensions:  # 仅保留与涡旋相关的变量
                    self.eddy_params[var_name] = var[:]

            # 遍历每个涡旋，初始化表层（depth_level=0）单元格
            for eddy_id in tqdm(range(num_eddies), desc="初始化表层涡旋单元格（depth_level=0）"):
                # 过滤条件1：涡旋有效半径<罗斯贝半径，跳过处理
                if f['effective_radius'][eddy_id] < Rossby_Radius:
                    continue

                # 步骤1：初始化涡旋中心的表层单元格（depth_level=0）
                center_lat = eddy_lat[eddy_id]
                center_lon = eddy_lon[eddy_id]
                center_hex = self._get_hex_code(center_lat, center_lon)  # 中心单元格哈希键
                self._init_hex_cell(center_hex, 0, center_lat, center_lon)  # 初始化表层单元格
                self.grid[(center_hex, 0)]['eddy_id'] = eddy_id  # 分配涡旋ID
                self.grid[(center_hex, 0)]['is_cyclonic'] = (eddy_type[eddy_id] == 1)  # 标记涡旋类型

                # 步骤2：初始化涡旋边界及边界邻域的表层单元格（depth_level=0）
                covered_hexes = {center_hex}  # 已覆盖的单元格（去重）
                for k in range(boundary_sample_count):
                    # 边界经纬度转换（原始数据可能放大100倍，需还原）
                    b_lat = boundary_lat[eddy_id, k] / 100
                    b_lon = boundary_lon[eddy_id, k] / 100

                    # 过滤无效经纬度（超出[-90,90]纬度或[-180,180]经度）
                    if not (-90 <= b_lat <= 90) or not (-180 <= b_lon <= 180):
                        continue

                    # 初始化边界单元格
                    b_hex = self._get_hex_code(b_lat, b_lon)
                    covered_hexes.add(b_hex)
                    self._init_hex_cell(b_hex, 0, b_lat, b_lon)
                    self.grid[(b_hex, 0)]['eddy_id'] = eddy_id
                    self.grid[(b_hex, 0)]['is_cyclonic'] = (eddy_type[eddy_id] == 1)

                    # 初始化边界单元格的1层邻域（保证表层覆盖完整）
                    for nbr_hex in self._get_neighbors(b_hex):
                        # 解码邻域单元格的经纬度
                        nbr_code = PyDGGS.DGGSElementCode()
                        nbr_code.row, nbr_code.col, nbr_code.basePartition = nbr_hex
                        nbr_code.level = self.hex_resolution
                        nbr_code.elementType = element_type
                        nbr_coord = PyDGGS.GeoCoord()
                        kernel.decode(nbr_code, nbr_coord)

                        # 初始化邻域单元格
                        self._init_hex_cell(nbr_hex, 0, nbr_coord.lat, nbr_coord.lon)
                        self.grid[(nbr_hex, 0)]['eddy_id'] = eddy_id
                        covered_hexes.add(nbr_hex)

                # 步骤3：计算表层（depth_level=0）涡旋轮廓、中心与半径
                contour_hexes, contour_center, radius = self._calculate_contour(
                    covered_hexes, depth_level=0, eddy_id=eddy_id
                )

                # 整理表层涡旋单元格信息（哈希键、纬度、经度）
                eddy_cells_info = [
                    (hex_code, self.grid[(hex_code, 0)]['lat'], self.grid[(hex_code, 0)]['lon'])
                    for hex_code in covered_hexes
                ]

                # 保存涡旋基础信息（含表层深度层级数据）
                self.eddy_cells[eddy_id] = {
                    'center_hex': center_hex,  # 表层中心单元格哈希键
                    'center_lat': center_lat,  # 表层中心纬度（原始数据）
                    'center_lon': center_lon,  # 表层中心经度（原始数据）
                    'cells': eddy_cells_info,  # 表层涡旋单元格列表
                    'is_cyclonic': (eddy_type[eddy_id] == 1),  # 涡旋类型（1=气旋）
                    'depth_info': {
                        0: {  # 表层（depth_level=0）详细信息
                            'center_hex': center_hex,
                            'center_lon': contour_center[0] if contour_center else center_lon,  # 计算的中心经度（优先）
                            'center_lat': contour_center[1] if contour_center else center_lat,  # 计算的中心纬度（优先）
                            'cells': eddy_cells_info,
                            'contour': contour_hexes,  # 表层轮廓单元格列表
                            'radius': radius  # 表层最大半径（单位：米）
                        }
                    }
                }

    def _verify_eddy_rotation(self, eddy_id, depth_level, hex_codes):
        """
        验证指定深度层级下，涡旋单元格集合的旋转特性（是否符合气旋/反气旋物理意义）

        核心验证条件（共5项）：
            1. 数据质量：有效速度点≥4个（否则无法判断旋转方向）
            2. 速度结构连续性：相邻点速度比≤1.5*Sv（避免速度突变）
            3. 旋转方向正确性：气旋→逆时针（角度差多为正），反气旋→顺时针（角度差多为负）
            4. 极端角度差过滤：角度差绝对值>1.5*Sa的点占比≤20%（避免方向突变）
            5. 方向分布均匀性：8个角度区间的点占比≤40%（避免单点聚集）
            6. 旋转对称性：旋转偏移后的角度误差≤π/2（保证对称结构）

        Args:
            eddy_id (int): 涡旋ID
            depth_level (int): 深度层级索引（0~29，直接作为速度场深度维度索引）
            hex_codes (list): 目标深度层级的涡旋单元格哈希键列表

        Returns:
            bool: 验证结果（True=通过，False=失败）
        """
        # 步骤1：提取目标深度层级的速度数据（东向u、北向v）
        u_values = []  # 东向速度列表
        v_values = []  # 北向速度列表
        for hex_code in hex_codes:
            # 获取单元格数据（若不存在则跳过）
            cell = self.grid.get((hex_code, depth_level), {})
            lat, lon = cell.get('lat'), cell.get('lon')
            if lat is None or lon is None:
                continue

            # 读取速度场的经纬度网格（用于最近邻插值）
            vel_lat = self.velocity_data['latitude'][:]
            vel_lon = self.velocity_data['longitude'][:]

            # 最近邻插值：找到速度场中最接近单元格经纬度的网格索引
            vel_i = int(np.abs(vel_lat - lat).argmin())  # 纬度索引
            vel_j = int(np.abs(vel_lon - lon).argmin())  # 经度索引

            try:
                # 读取原始速度数据（速度场维度：time, depth, lat, lon，time固定取第0维）
                u_raw = self.velocity_data['uo'][0, depth_level, vel_i, vel_j]  # 东向原始速度
                v_raw = self.velocity_data['vo'][0, depth_level, vel_i, vel_j]  # 北向原始速度

                # 过滤掩码数据（掩码表示无效值）
                if np.ma.is_masked(u_raw) or np.ma.is_masked(v_raw):
                    continue

                # 应用速度场的缩放因子和偏移量（还原真实速度值）
                u = u_raw * self.velocity_data['uo'].attrs.get('scale_factor', 1.0) + \
                    self.velocity_data['uo'].attrs.get('add_offset', 0.0)
                v = v_raw * self.velocity_data['vo'].attrs.get('scale_factor', 1.0) + \
                    self.velocity_data['vo'].attrs.get('add_offset', 0.0)

                u_values.append(u)
                v_values.append(v)
            except Exception:
                # 捕获索引越界、变量不存在等异常，跳过该单元格
                continue

        # 验证条件1：数据质量（有效速度点≥4个，否则无法判断旋转方向）
        velocities = np.column_stack((u_values, v_values)) if len(u_values) >= 4 else np.array([])
        if len(velocities) < 4:
            self.failure_stats['data_quality'] += 1
            return False

        # 预处理：计算速度大小（模长）和方向角（弧度）
        vel_mags = np.linalg.norm(velocities, axis=1)  # 速度大小（每个点1个值）
        dir_angles = np.arctan2(velocities[:, 1], velocities[:, 0])  # 速度方向角（弧度）
        is_cyclonic = self.eddy_cells[eddy_id]['is_cyclonic']  # 当前涡旋类型（气旋/反气旋）
        n_points = len(velocities)  # 有效速度点数量

        # 验证条件2：速度结构连续性（相邻点速度比≤1.5*Sv）
        max_vel_ratio = self.params['Sv'] * 1.5  # 最大允许速度比
        valid_ratio_count = 0  # 满足速度比条件的点数
        for i in range(n_points):
            next_idx = (i + 1) % n_points  # 下一个点的索引（循环）
            # 避免除以0（速度大小接近0时跳过）
            if vel_mags[next_idx] < 1e-6:
                continue
            # 计算当前点与下一点的速度比
            ratio = vel_mags[i] / vel_mags[next_idx]
            if (1 / max_vel_ratio) <= ratio <= max_vel_ratio:
                valid_ratio_count += 1
        # 满足条件的点数占比<80%，验证失败
        if valid_ratio_count / n_points < 0.8:
            self.failure_stats['velocity_structure'] += 1
            return False

        # 验证条件3：旋转方向正确性（气旋逆时针/反气旋顺时针）
        sa_rad = (self.params['Sa'] * np.pi / 180) * 1.5  # 最大允许角度差（弧度，放大1.5倍）
        sae_rad = (self.params['Sae'] * np.pi / 180) * 2.0  # 方向容忍阈值（弧度，放大2倍）

        # 计算相邻点的方向角差（归一化到[-π, π]）
        angle_diffs = []
        for i in range(n_points):
            next_idx = (i + 1) % n_points
            diff = dir_angles[next_idx] - dir_angles[i]
            # 角度差归一化：确保在[-π, π]范围内
            diff = (diff + np.pi) % (2 * np.pi) - np.pi
            angle_diffs.append(diff)
        angle_diffs = np.array(angle_diffs)
        mean_diff = np.mean(angle_diffs)  # 角度差平均值

        # 气旋：角度差> -sae_rad的点占比≥60% 且 平均值≥-0.1（逆时针趋势）
        if is_cyclonic:
            valid_diffs = angle_diffs[angle_diffs > -sae_rad]
            if len(valid_diffs) / n_points < 0.6 or mean_diff < -0.1:
                self.failure_stats['rotation_direction'] += 1
                return False
        # 反气旋：角度差< sae_rad的点占比≥60% 且 平均值≤0.1（顺时针趋势）
        else:
            valid_diffs = angle_diffs[angle_diffs < sae_rad]
            if len(valid_diffs) / n_points < 0.6 or mean_diff > 0.1:
                self.failure_stats['rotation_direction'] += 1
                return False

        # 验证条件4：极端角度差过滤（角度差绝对值>sa_rad的点占比≤20%）
        extreme_diffs = angle_diffs[np.abs(angle_diffs) > sa_rad]
        if len(extreme_diffs) / n_points > 0.2:
            self.failure_stats['symmetry'] += 1
            return False

        # 验证条件5：方向分布均匀性（8个角度区间的最大占比≤40%）
        # 将方向角[-π, π]分为8个区间，统计每个区间的点数
        angle_hist, _ = np.histogram(dir_angles, bins=8, range=(-np.pi, np.pi))
        # 最大区间点数占比>40%，说明方向分布不均匀
        if np.max(angle_hist) / n_points > 0.4:
            self.failure_stats['symmetry'] += 1
            return False

        # 验证条件6：旋转对称性（旋转偏移后的角度误差≤π/2）
        rotation_sim = []  # 旋转对称性误差列表
        for offset in range(1, n_points // 2):
            # 将方向角列表滚动offset个位置（模拟旋转）
            shifted_angles = np.roll(dir_angles, offset)
            # 计算旋转后的角度误差（取最小角度差，避免跨0点误差）
            angle_errors = np.minimum(
                np.abs(dir_angles - shifted_angles),
                2 * np.pi - np.abs(dir_angles - shifted_angles)
            )
            rotation_sim.append(np.mean(angle_errors))  # 平均误差

        # 最小平均误差>π/2，说明旋转对称性差
        if np.min(rotation_sim) > np.pi / 2:
            self.failure_stats['symmetry'] += 1
            return False

        # 综合评分（需≥0.5才通过验证）
        angle_score = len(valid_diffs) / n_points  # 旋转方向评分（0~1）
        # 方向分布熵（0~log2(8)=3），归一化到0~1
        entropy = -np.sum([(c / n_points) * np.log2(c / n_points + 1e-10) for c in angle_hist if c > 0])
        deadzone_score = entropy / np.log2(8)  # 方向分布均匀性评分（0~1）
        # 旋转对称性评分（0~1，误差越小评分越高）
        symmetry_score = 1 - (np.min(rotation_sim) / (np.pi / 2))
        # 总评分（权重：旋转方向50%，分布均匀性30%，对称性20%）
        total_score = angle_score * 0.5 + deadzone_score * 0.3 + symmetry_score * 0.2

        return total_score >= 0.5

    def fuse_velocity_data(self, velocity_file, max_eddies=None):
        """
        融合三维速度场数据，实现涡旋的分层外扩提取与验证

        核心逻辑（业务需求）：
            1. 深度层级（0~29）与速度场深度维度一一对应，无额外映射
            2. 从表层（depth_level=0）开始，逐层向下外扩提取涡旋
            3. 每层外扩前需通过旋转特性验证，失败则停在上一层
            4. 下层涡旋轮廓半径≤上层半径（避免下层过度扩张）
            5. 提取完成后，过滤虚假涡旋（深度层级总数<3）

        Args:
            velocity_file (str): 三维速度场数据NetCDF文件路径
            max_eddies (int, optional): 最大处理涡旋数（用于测试，None=处理全部）
        """
        # 读取速度场NetCDF文件
        with nc.Dataset(velocity_file, 'r') as f:
            self.velocity_data = f  # 保存速度场数据对象

            # 筛选待处理的涡旋ID（支持限制最大数量，用于快速测试）
            eddy_ids = list(self.eddy_cells.keys())
            if max_eddies is not None and max_eddies > 0:
                eddy_ids = eddy_ids[:max_eddies]
                print(f"启用快速测试模式：仅处理前 {max_eddies} 个涡旋")

            # 遍历每个涡旋，处理分层外扩（从depth_level=1开始，表层已初始化）
            for eddy_id in tqdm(eddy_ids, desc="涡旋分层外扩与速度场融合"):
                # 跳过已标记为虚假涡旋的ID
                if eddy_id in self.false_eddies:
                    continue

                # 逐层处理深度层级（从1到29，表层0已处理）
                for depth_level in range(1, self.depth_levels):
                    prev_depth_level = depth_level - 1  # 上一层深度层级（用于约束当前层）

                    # 业务需求3：上一层无数据时，停止当前涡旋的深度扩展
                    if prev_depth_level not in self.eddy_cells[eddy_id]['depth_info']:
                        break

                    # 步骤1：获取上一层信息（作为当前层的约束）
                    prev_info = self.eddy_cells[eddy_id]['depth_info'][prev_depth_level]
                    prev_center_hex = prev_info['center_hex']  # 上一层中心单元格
                    upper_radius = prev_info['radius']  # 上一层半径（当前层≤该值）

                    # 步骤2：初始外扩（中心+1层邻域），第一次验证
                    initial_valid = {prev_center_hex}  # 初始有效单元格（上一层中心）
                    # 初始外围：上一层中心的1层邻域
                    initial_perimeter = self._get_neighbors(prev_center_hex, radius=1)
                    initial_perimeter = list(set(initial_perimeter) - initial_valid)  # 去重（排除中心）
                    initial_candidates = list(initial_valid.union(initial_perimeter))  # 初始候选单元格

                    # 第一次验证失败：当前层无数据，停止扩展
                    if not self._verify_eddy_rotation(eddy_id, depth_level, initial_candidates):
                        self.failure_stats['rotation_direction'] += 1
                        break

                    # 步骤3：分层外扩循环（业务需求2：层层外扩；业务需求5：半径≤上层）
                    current_valid = initial_valid  # 当前层已验证通过的单元格
                    current_perimeter = initial_perimeter  # 当前层待验证的外围单元格
                    layer_continue = True  # 是否继续外扩的标志

                    while current_perimeter and layer_continue:
                        # a. 生成新候选单元格（外围单元格的1层邻域，排除已有效单元格）
                        new_candidates = set()
                        for hex_code in current_perimeter:
                            new_candidates.update(self._get_neighbors(hex_code, radius=1))
                        new_candidates = new_candidates - current_valid  # 去重（排除已有效）
                        new_candidates = list(new_candidates)
                        if not new_candidates:
                            break  # 无新候选，停止外扩

                        # b. 验证新候选的旋转特性（需结合已有有效单元格）
                        all_candidates = list(current_valid.union(new_candidates))
                        if not self._verify_eddy_rotation(eddy_id, depth_level, all_candidates):
                            break  # 验证失败，停止外扩

                        # c. 业务需求5：当前层半径≤上层半径，否则强制停止
                        temp_hexes = list(current_valid.union(new_candidates))
                        _, _, current_radius = self._calculate_contour(temp_hexes, depth_level, eddy_id)
                        if current_radius > upper_radius:
                            layer_continue = False
                            break

                        # d. 更新有效单元格和外围（外扩一步）
                        current_valid.update(new_candidates)  # 新候选加入有效集合
                        # 新外围：新候选的1层邻域，排除已有效单元格
                        new_perimeter = set()
                        for hex_code in new_candidates:
                            new_perimeter.update(self._get_neighbors(hex_code, radius=1))
                        current_perimeter = list(new_perimeter - current_valid)

                    # 步骤4：初始化当前层有效单元格，并填充速度数据
                    current_cells = []  # 当前层涡旋单元格列表（哈希键、纬度、经度）
                    for hex_code in current_valid:
                        # 解码单元格经纬度
                        code = PyDGGS.DGGSElementCode()
                        code.row, code.col, code.basePartition = hex_code
                        code.level = self.hex_resolution
                        code.elementType = element_type
                        coord = PyDGGS.GeoCoord()
                        kernel.decode(code, coord)
                        lat = coord.lat
                        lon = coord.lon

                        # 初始化当前层单元格
                        self._init_hex_cell(hex_code, depth_level, lat, lon)
                        current_cells.append((hex_code, lat, lon))

                        # 填充速度数据（最近邻插值）
                        vel_lat = self.velocity_data['latitude'][:]
                        vel_lon = self.velocity_data['longitude'][:]
                        vel_i = int(np.abs(vel_lat - lat).argmin())  # 纬度索引
                        vel_j = int(np.abs(vel_lon - lon).argmin())  # 经度索引

                        try:
                            # 读取原始速度数据（time=0，depth=当前depth_level）
                            u_raw = self.velocity_data['uo'][0, depth_level, vel_i, vel_j]
                            v_raw = self.velocity_data['vo'][0, depth_level, vel_i, vel_j]

                            # 过滤掩码数据，填充真实速度值
                            if not np.ma.is_masked(u_raw) and not np.ma.is_masked(v_raw):
                                u = u_raw * self.velocity_data['uo'].attrs.get('scale_factor', 1.0) + \
                                    self.velocity_data['uo'].attrs.get('add_offset', 0.0)
                                v = v_raw * self.velocity_data['vo'].attrs.get('scale_factor', 1.0) + \
                                    self.velocity_data['vo'].attrs.get('add_offset', 0.0)
                                # 更新单元格速度和涡旋ID
                                self.grid[(hex_code, depth_level)]['u'] = u
                                self.grid[(hex_code, depth_level)]['v'] = v
                                self.grid[(hex_code, depth_level)]['eddy_id'] = eddy_id
                        except Exception:
                            # 捕获异常，跳过该单元格的速度填充（保留NaN）
                            continue

                    # 步骤5：计算当前层轮廓、中心与半径，保存深度信息
                    contour_hexes, center, final_radius = self._calculate_contour(
                        [hex_code for hex_code, _, _ in current_cells], depth_level, eddy_id
                    )
                    # 中心计算失败时，沿用前一层中心
                    if not center:
                        center = (prev_info['center_lon'], prev_info['center_lat'])
                    # 计算当前层中心单元格哈希键
                    current_center_hex = self._get_hex_code(center[1], center[0])

                    # 保存当前层深度信息
                    self.eddy_cells[eddy_id]['depth_info'][depth_level] = {
                        'center_hex': current_center_hex,
                        'center_lon': center[0],
                        'center_lat': center[1],
                        'cells': current_cells,
                        'contour': contour_hexes,
                        'radius': final_radius
                    }

                # 步骤6：业务需求4：过滤虚假涡旋（深度层级总数<3）
                # 深度层级从0开始，总层数=最大深度层级+1（如0/1/2 → 3层）
                depth_keys = self.eddy_cells[eddy_id]['depth_info'].keys()
                max_depth_level = max(depth_keys) if depth_keys else -1
                total_layers = max_depth_level + 1
                if total_layers < 3:
                    self.false_eddies.append(eddy_id)

            # 输出处理统计信息
            print(f"\n=== 涡旋分层处理统计 ===")
            print(f"- 总处理涡旋数：{len(eddy_ids)}")
            print(f"- 虚假涡旋数（深度层级<3）：{len(self.false_eddies)}")
            print(f"- 验证失败原因统计：{self.failure_stats}")

    def _haversine_distance(self, lat1, lon1, lat2, lon2):
        """
        哈维正弦公式（Haversine Formula）：计算地球表面两点间的最短距离（球面距离）

        Args:
            lat1 (float): 第一点纬度（单位：度）
            lon1 (float): 第一点经度（单位：度）
            lat2 (float): 第二点纬度（单位：度）
            lon2 (float): 第二点经度（单位：度）

        Returns:
            float: 两点间距离（单位：米）
        """
        R = 6371000  # 地球平均半径（单位：米）
        # 转换角度为弧度（三角函数要求）
        phi1, phi2 = math.radians(lat1), math.radians(lat2)
        dphi = math.radians(lat2 - lat1)  # 纬度差（弧度）
        dlambda = math.radians(lon2 - lon1)  # 经度差（弧度）

        # 哈维正弦公式计算
        a = math.sin(dphi / 2) ** 2 + math.cos(phi1) * math.cos(phi2) * math.sin(dlambda / 2) ** 2
        c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))  # 圆心角（弧度）

        return R * c  # 距离=半径×圆心角

    def save_results(self, output_dir, base_name, eddy_file_path):
        """
        保存涡旋提取结果（二维表层数据+三维全深度数据）

        输出规则：
            1. 按涡旋类型分类目录：气旋→2DC（二维）/3DC（三维），反气旋→2DA/3DA
            2. 二维文件：仅保存表层（depth_level=0）数据，含原始涡旋元数据
            3. 三维文件：保存所有有效深度层级数据，含速度、经纬度、轮廓信息

        Args:
            output_dir (str): 结果输出根目录
            base_name (str): 输出文件基础名（从涡旋文件名提取，不含后缀）
            eddy_file_path (str): 涡旋原始数据文件路径（用于复制元数据）
        """
        # 步骤1：创建输出目录（按涡旋类型分类）
        two_d_folder = os.path.join(output_dir, "2DC" if self.is_cyclonic else "2DA")  # 二维目录
        three_d_folder = os.path.join(output_dir, "3DC" if self.is_cyclonic else "3DA")  # 三维目录
        os.makedirs(two_d_folder, exist_ok=True)  # 不存在则创建，存在则忽略
        os.makedirs(three_d_folder, exist_ok=True)

        # 步骤2：筛选有效涡旋（排除虚假涡旋）
        valid_eddies = [eid for eid in self.eddy_cells if eid not in self.false_eddies]
        num_valid = len(valid_eddies)
        if num_valid == 0:
            print(f"警告：无有效涡旋数据（{base_name}），跳过保存")
            return

        # 步骤3：读取原始涡旋文件的元数据（用于复制到输出文件）
        with nc.Dataset(eddy_file_path, 'r') as orig_file:
            nb_sample = len(orig_file.dimensions['NbSample'])  # 边界采样点数
            global_attrs = dict(orig_file.attrs)  # 全局属性（如数据来源、创建时间）
            var_attrs = {var: dict(orig_file[var].attrs) for var in orig_file.variables}  # 变量属性

        # ---------------------- 保存二维数据（仅表层：depth_level=0）----------------------
        two_d_output = os.path.join(two_d_folder, f"{base_name}_2D.nc")
        with nc.Dataset(two_d_output, 'w') as f:
            # 定义维度（与原始涡旋文件一致）
            f.dimensions['obs'] = num_valid  # 有效涡旋数
            f.dimensions['NbSample'] = nb_sample  # 边界采样点数

            # 复制全局属性（保持数据溯源性）
            for attr_name, attr_val in global_attrs.items():
                f.attrs[attr_name] = attr_val

            # 定义并填充二维变量（原始涡旋元数据）
            # 注：变量列表需与原始涡旋文件的变量名匹配，可根据实际数据调整
            two_d_vars = [
                ('amplitude', ('obs',), np.uint16),
                ('effective_area', ('obs',), np.float32),
                ('effective_contour_height', ('obs',), np.float32),
                ('effective_contour_latitude', ('obs', 'NbSample'), np.int16),
                ('effective_contour_longitude', ('obs', 'NbSample'), np.int16),
                ('effective_contour_shape_error', ('obs',), np.uint8),
                ('effective_radius', ('obs',), np.uint16),
                ('inner_contour_height', ('obs',), np.float32),
                ('latitude', ('obs',), np.float32),
                ('latitude_max', ('obs',), np.float32),
                ('longitude', ('obs',), np.float32),
                ('longitude_max', ('obs',), np.float32),
                ('num_contours', ('obs',), np.uint16),
                ('num_point_e', ('obs',), np.uint16),
                ('num_point_s', ('obs',), np.uint16),
                ('speed_area', ('obs',), np.float32),
                ('speed_average', ('obs',), np.uint16),
                ('speed_contour_height', ('obs',), np.float32),
                ('speed_contour_latitude', ('obs', 'NbSample'), np.int16),
                ('speed_contour_longitude', ('obs', 'NbSample'), np.int16),
                ('speed_contour_shape_error', ('obs',), np.uint8),
                ('speed_radius', ('obs',), np.uint16),
                ('time', ('obs',), np.uint32),
                ('uavg_profile', ('obs', 'NbSample'), np.uint16),
            ]

            # 从原始文件读取数据，填充到二维输出文件
            with nc.Dataset(eddy_file_path, 'r') as orig_file:
                for name, dims, dtype in two_d_vars:
                    # 创建变量并复制原始属性
                    var = f.createVariable(name, dtype, dims)
                    for attr_name, attr_val in var_attrs[name].items():
                        var.attrs[attr_name] = attr_val
                    # 按有效涡旋ID映射数据（新索引→原始索引）
                    for new_idx, old_idx in enumerate(valid_eddies):
                        var[new_idx] = orig_file[name][old_idx]

        print(f"二维表层数据已保存：{two_d_output}")

        # ---------------------- 保存三维数据（所有有效深度层级）----------------------
        three_d_output = os.path.join(three_d_folder, f"{base_name}_3D.nc")
        with nc.Dataset(three_d_output, 'w') as f:
            # 定义维度（新增depth_level维度，与速度场一致）
            f.dimensions['obs'] = num_valid  # 有效涡旋数
            f.dimensions['NbSample'] = nb_sample  # 边界采样点数
            f.dimensions['depth_level'] = self.depth_levels  # 深度层级（0~29）
            f.dimensions['max_points'] = 1000  # 每层最大单元格数（可根据数据调整）
            f.dimensions['contour_points'] = 16  # 每层轮廓点数（与_calculate_contour一致）

            # 复制全局属性
            for attr_name, attr_val in global_attrs.items():
                f.attrs[attr_name] = attr_val

            # 步骤1：复制二维变量（与二维文件一致，保持元数据完整性）
            with nc.Dataset(eddy_file_path, 'r') as orig_file:
                for name, dims, dtype in two_d_vars:
                    var = f.createVariable(name, dtype, dims)
                    for attr_name, attr_val in var_attrs[name].items():
                        var.attrs[attr_name] = attr_val
                    for new_idx, old_idx in enumerate(valid_eddies):
                        var[new_idx] = orig_file[name][old_idx]

            # 步骤2：定义三维特有变量（按深度层级组织的数据）
            three_d_vars = [
                # 每层中心与半径
                ('center_latitude', ('obs', 'depth_level'), np.float32),  # 中心纬度
                ('center_longitude', ('obs', 'depth_level'), np.float32),  # 中心经度
                ('radius', ('obs', 'depth_level'), np.float32),  # 最大半径（米）
                # 每层单元格的速度与经纬度
                ('u_velocity', ('obs', 'depth_level', 'max_points'), np.float32),  # 东向速度
                ('v_velocity', ('obs', 'depth_level', 'max_points'), np.float32),  # 北向速度
                ('point_latitude', ('obs', 'depth_level', 'max_points'), np.float32),  # 单元格纬度
                ('point_longitude', ('obs', 'depth_level', 'max_points'), np.float32),  # 单元格经度
                # 每层轮廓的经纬度
                ('contour_latitude', ('obs', 'depth_level', 'contour_points'), np.float32),  # 轮廓纬度
                ('contour_longitude', ('obs', 'depth_level', 'contour_points'), np.float32),  # 轮廓经度
            ]

            # 创建三维变量并初始化为NaN（无效值用NaN表示）
            for name, dims, dtype in three_d_vars:
                var = f.createVariable(name, dtype, dims)
                var[:] = np.nan  # 初始填充NaN

            # 步骤3：填充三维数据（按涡旋→深度层级的顺序）
            for e_idx, eddy_id in enumerate(valid_eddies):
                depth_info = self.eddy_cells[eddy_id]['depth_info']  # 当前涡旋的所有深度信息

                # 遍历所有深度层级（0~29）
                for depth_level in range(self.depth_levels):
                    # 跳过无数据的深度层级（保留NaN）
                    if depth_level not in depth_info:
                        continue

                    info = depth_info[depth_level]  # 当前深度层级的详细信息

                    # 填充中心与半径
                    f['center_latitude'][e_idx, depth_level] = info['center_lat']
                    f['center_longitude'][e_idx, depth_level] = info['center_lon']
                    f['radius'][e_idx, depth_level] = info['radius']

                    # 填充单元格数据（经纬度+速度）：限制最大数量为max_points
                    cells = info['cells'][:f.dimensions['max_points'].size]
                    for p_idx, (hex_code, lat, lon) in enumerate(cells):
                        # 单元格经纬度
                        f['point_latitude'][e_idx, depth_level, p_idx] = lat
                        f['point_longitude'][e_idx, depth_level, p_idx] = lon
                        # 单元格速度（从grid中读取，无则保留NaN）
                        cell_data = self.grid.get((hex_code, depth_level), {})
                        f['u_velocity'][e_idx, depth_level, p_idx] = cell_data.get('u', np.nan)
                        f['v_velocity'][e_idx, depth_level, p_idx] = cell_data.get('v', np.nan)

                    # 填充轮廓数据（经纬度）：限制数量为contour_points，不足则用NaN填充
                    contour_hexes = info['contour'][:f.dimensions['contour_points'].size]
                    contour_points = []  # 存储轮廓点（lat, lon）

                    # 解码轮廓单元格的经纬度
                    for hex_code in contour_hexes:
                        code = PyDGGS.DGGSElementCode()
                        code.row, code.col, code.basePartition = hex_code
                        code.level = self.hex_resolution
                        code.elementType = element_type
                        coord = PyDGGS.GeoCoord()
                        kernel.decode(code, coord)
                        contour_points.append((coord.lat, coord.lon))

                    # 不足16个轮廓点时，用NaN填充
                    while len(contour_points) < f.dimensions['contour_points'].size:
                        contour_points.append((np.nan, np.nan))

                    # 写入轮廓数据
                    for c_idx, (lat, lon) in enumerate(contour_points):
                        f['contour_latitude'][e_idx, depth_level, c_idx] = lat
                        f['contour_longitude'][e_idx, depth_level, c_idx] = lon

        print(f"三维全深度数据已保存：{three_d_output}")


def extract_date_identifier(filename):
    """
    从文件名中提取日期标识（8位数字，如20240101）

    Args:
        filename (str): 文件名（含或不含路径）

    Returns:
        str or None: 提取的8位日期字符串，无匹配则返回None
    """
    # 正则表达式匹配8位数字（日期）
    match = re.search(r'\d{8}', filename)
    return match.group() if match else None


def find_matching_file(directory, identifier):
    """
    在指定目录中查找含指定标识（如日期）的NetCDF文件（用于匹配涡旋与速度场文件）

    Args:
        directory (str): 查找目录路径
        identifier (str): 匹配标识（如8位日期字符串）

    Returns:
        str or None: 匹配文件的完整路径，无匹配则返回None
    """
    # 输入无效时返回None
    if not identifier or not os.path.exists(directory):
        return None

    # 遍历目录下的文件，匹配含标识且后缀为.nc的文件
    for filename in os.listdir(directory):
        if identifier in filename and filename.endswith('.nc'):
            return os.path.join(directory, filename)

    return None


def process_eddy_files(eddy_dir, velocity_dir, output_dir, max_eddies_per_file=None):
    """
    批量处理涡旋文件：自动匹配速度场文件→处理涡旋→保存结果

    处理流程：
        1. 遍历涡旋目录下的所有NetCDF文件
        2. 从涡旋文件名提取日期标识，匹配对应的速度场文件
        3. 根据涡旋文件名判断涡旋类型（气旋/反气旋）
        4. 初始化EddyHybridIntegration对象，完成表层初始化与速度场融合
        5. 保存二维（表层）与三维（全深度）结果

    Args:
        eddy_dir (str): 涡旋原始数据目录路径
        velocity_dir (str): 三维速度场数据目录路径
        output_dir (str): 结果输出根目录路径
        max_eddies_per_file (int, optional): 每文件最大处理涡旋数（测试用，None=全部）
    """
    # 创建输出根目录（不存在则创建）
    os.makedirs(output_dir, exist_ok=True)

    # 筛选涡旋目录下的NetCDF文件（仅处理.nc后缀）
    eddy_files = [f for f in os.listdir(eddy_dir) if f.endswith('.nc')]
    print(f"在涡旋目录 {eddy_dir} 中发现 {len(eddy_files)} 个NetCDF文件")

    # 遍历每个涡旋文件，批量处理
    for eddy_file in eddy_files:
        eddy_path = os.path.join(eddy_dir, eddy_file)
        print(f"\n=== 开始处理涡旋文件：{eddy_file} ===")

        # 步骤1：提取日期标识（用于匹配速度场文件）
        date_id = extract_date_identifier(eddy_file)
        if not date_id:
            print(f"警告：无法从文件名 {eddy_file} 提取8位日期标识，跳过该文件")
            continue

        # 步骤2：匹配对应的速度场文件（按日期标识）
        velocity_file = find_matching_file(velocity_dir, date_id)
        if not velocity_file:
            print(f"警告：未找到与日期 {date_id} 匹配的速度场文件，跳过该文件")
            continue
        print(f"成功匹配速度场文件：{os.path.basename(velocity_file)}")

        # 步骤3：判断涡旋类型（文件名含"anticyclonic"则为反气旋，否则为气旋）
        is_cyclonic = not ("anticyclonic" in eddy_file.lower())
        print(f"涡旋类型判定：{'气旋（Cyclonic）' if is_cyclonic else '反气旋（Anticyclonic）'}")

        # 步骤4：初始化涡旋处理器并处理数据
        # 注：hex_resolution默认与全局dggs_level一致，可根据需求调整
        integrator = EddyHybridIntegration(is_cyclonic=is_cyclonic, hex_resolution=dggs_level)
        # 初始化表层（depth_level=0）涡旋单元格与ID分配
        integrator.assign_eddy_ids(eddy_path, is_cyclonic)
        # 融合速度场数据，分层外扩与验证
        integrator.fuse_velocity_data(velocity_file, max_eddies=max_eddies_per_file)

        # 步骤5：保存处理结果（二维+三维）
        # 输出文件基础名：去掉涡旋文件名的.nc后缀
        base_name = os.path.basename(eddy_file).replace('.nc', '')
        integrator.save_results(output_dir, base_name, eddy_path)

    print(f"\n=== 所有涡旋文件处理完成 ===")


if __name__ == "__main__":
    """
    主程序入口：配置路径并启动批量处理

    路径配置说明：
        - eddy_dir: 涡旋原始数据目录（存放涡旋NetCDF文件）
        - velocity_dir: 三维速度场数据目录（存放速度场NetCDF文件）
        - output_dir: 结果输出目录（自动创建2DC/2DA/3DC/3DA子目录）
        - max_eddies_per_file: 每文件最大处理涡旋数（测试用，None=处理全部）
    """
    # ---------------------- 路径配置（请根据实际数据路径修改）----------------------
    eddy_dir = "./EddyExtract/Ay"  # 涡旋原始数据目录
    velocity_dir = "./3DV"  # 三维速度场数据目录
    output_dir = "./fused_results_hex"  # 结果输出根目录
    max_eddies_per_file = None  # 测试时可设为10（仅处理前10个涡旋）

    # ---------------------- 启动批量处理 ----------------------
    process_eddy_files(
        eddy_dir=eddy_dir,
        velocity_dir=velocity_dir,
        output_dir=output_dir,
        max_eddies_per_file=max_eddies_per_file
    )